# -*- coding: utf-8 -*-
"""
Created on Wed May  2 17:06:25 2018
Tests 3 phase power flow algorithm
@author: sghosh
"""
import os

import numpy as np
import pytest
import copy

from pandapower.networks import create_cigre_network_mv
from pandapower import pp_dir
from pandapower import replace_line_by_impedance
from pandapower.auxiliary import get_free_id
from pandapower.create import create_empty_network, create_bus, create_ext_grid, create_line, create_asymmetric_load, \
    create_line_from_parameters, create_load, create_transformer_from_parameters, create_sgen, create_asymmetric_sgen, \
    create_shunt, create_gen, create_switch
from pandapower.file_io import from_json
from pandapower.pf.runpp_3ph import runpp_3ph
from pandapower.run import runpp
from pandapower.std_types import create_std_type, add_zero_impedance_parameters
from pandapower.test.consistency_checks import runpp_3ph_with_consistency_checks, runpp_with_consistency_checks, \
    trafo_currents_consistent_3ph
from pandapower.test.loadflow.PF_Results import get_PF_Results
from pandapower.toolbox import dataframes_equal
from pandapower.test.conftest import result_test_network
from pandapower.test.helper_functions import add_grid_connection
from pandapower.test.loadflow.test_runpp import get_isolated


@pytest.fixture
def net():
    v_base = 110  # 110kV Base Voltage
    k_va_base = 100  # 100 MVA
    net = create_empty_network(sn_mva=k_va_base)
    create_bus(net, vn_kv=v_base, index=1)
    create_bus(net, vn_kv=v_base, index=5)
    create_ext_grid(net, bus=1, vm_pu=1.0, s_sc_max_mva=5000, rx_max=0.1,
                    r0x0_max=0.1, x0x_max=1.0)
    create_std_type(net, {"r0_ohm_per_km": 0.0848, "x0_ohm_per_km": 0.4649556,
                          "c0_nf_per_km": 230.6, "g0_us_per_km": 0, "max_i_ka": 0.963, "r_ohm_per_km": 0.0212,
                          "x_ohm_per_km": 0.1162389, "c_nf_per_km": 230},
                    "example_type")
    create_line(net, from_bus=1, to_bus=5, length_km=50.0, std_type="example_type")

    create_asymmetric_load(net, 5, p_a_mw=50, q_a_mvar=50, p_b_mw=10, q_b_mvar=15,
                           p_c_mw=10, q_c_mvar=5)
    return net


def check_it(net):
    bus_pp = np.abs(net.res_bus_3ph[['vm_a_pu', 'vm_b_pu', 'vm_c_pu']]
                    [~np.isnan(net.res_bus_3ph.vm_a_pu)].values)
    bus_pf = np.abs(np.array([[0.96742893, 1.01302766, 1.019784],
                              [0.74957533, 1.09137945, 1.05124282]]))
    assert np.max(np.abs(bus_pp - bus_pf)) < 1.1e-6
    line_pp = np.abs(net.res_line_3ph[~np.isnan(net.res_line_3ph.i_a_from_ka)]
                     [['i_a_from_ka', 'i_a_to_ka', 'i_b_from_ka', 'i_b_to_ka',
                       'i_c_from_ka', 'i_c_to_ka',
                       'p_a_from_mw', 'p_a_to_mw', 'q_a_from_mvar', 'q_a_to_mvar',
                       'p_b_from_mw', 'p_b_to_mw', 'q_b_from_mvar', 'q_b_to_mvar',
                       'p_c_from_mw', 'p_c_to_mw', 'q_c_from_mvar', 'q_c_to_mvar',
                       'loading_a_percent', 'loading_b_percent', 'loading_c_percent',
                       'loading_percent'
                       ]].values)
    line_pf = np.abs(np.array(
        [[1.34212045, 1.48537916, 0.13715552, 0.26009611,
          0.22838401, 0.1674634,
          55.70772301, (-49.999992954), 60.797262682, (-49.999959283),
          8.7799379802, (-9.9999996625), (-0.88093549983), (-15.000000238),
          9.3739293122, (-10.000000161), (-11.441663679), (-4.9999997418),
          154.2452, 27.00894, 23.71589,
          154.2452]]))
    assert np.max(np.abs(line_pp - line_pf)) < 1.1e-4


def test_2bus_network(net):
    # -o---o
    add_zero_impedance_parameters(net)
    runpp_3ph_with_consistency_checks(net)
    assert net['converged']
    check_it(net)


def test_2bus_network_single_isolated_busses(net):
    # -o---o o x
    create_bus(net, vn_kv=110)
    create_bus(net, vn_kv=110, in_service=False)
    add_zero_impedance_parameters(net)
    runpp_3ph_with_consistency_checks(net)
    assert net['converged']
    check_it(net)


def test_2bus_network_isolated_net_part(net):
    # -o---o o---o
    b1 = create_bus(net, vn_kv=110)
    b2 = create_bus(net, vn_kv=110)
    create_line(net, from_bus=b1, to_bus=b2, length_km=50.0, std_type="example_type")
    create_asymmetric_load(net, b2, p_a_mw=50, q_a_mvar=50, p_b_mw=10, q_b_mvar=15,
                           p_c_mw=10, q_c_mvar=5)
    add_zero_impedance_parameters(net)
    runpp_3ph_with_consistency_checks(net)
    assert net['converged']
    check_it(net)


def test_2bus_network_singel_oos_bus(net):
    # -o---x---o
    b1 = create_bus(net, vn_kv=110)
    net.bus.loc[5, "in_service"] = False
    create_line(net, from_bus=5, to_bus=b1, length_km=10.0, std_type="example_type")
    create_asymmetric_load(net, b1, p_a_mw=-5, q_a_mvar=5, p_b_mw=-1, q_b_mvar=1.5,
                           p_c_mw=-1, q_c_mvar=.5)
    add_zero_impedance_parameters(net)
    runpp_3ph_with_consistency_checks(net)
    assert net['converged']


def test_out_serv_load(net):
    # <-x--o------o
    add_zero_impedance_parameters(net)
    runpp_3ph_with_consistency_checks(net)
    assert net['converged']
    check_it(net)
    create_asymmetric_load(net, 5, p_a_mw=50, q_a_mvar=100, p_b_mw=29, q_b_mvar=38,
                           p_c_mw=10, q_c_mvar=5, in_service=False)
    runpp_3ph_with_consistency_checks(net)
    assert net['converged']
    check_it(net)


def test_2bus_network_one_of_two_ext_grids_oos(net):
    """ Out of service ext_grids should not affect the unbalanced powerflow"""
    create_ext_grid(net, bus=1, vm_pu=1.0, s_sc_max_mva=5000, rx_max=0.1,
                    r0x0_max=0.1, x0x_max=1.0, in_service=False)
    runpp_3ph(net)
    assert net['converged']


@pytest.mark.parametrize("init", ["auto", "results", "flat", "dc"])
@pytest.mark.parametrize("recycle", [None, {"bus_pq": True, "Ybus": True}, {"bus_pq": True, "Ybus": False}])
def test_4bus_network(init, recycle):
    v_base = 110  # 110kV Base Voltage
    mva_base = 100  # 100 MVA
    net = create_empty_network(sn_mva=mva_base)
    # =============================================================================
    # Main Program
    # =============================================================================
    busn = create_bus(net, vn_kv=v_base, name="busn")
    busk = create_bus(net, vn_kv=v_base, name="busk")
    busm = create_bus(net, vn_kv=v_base, name="busm")
    busp = create_bus(net, vn_kv=v_base, name="busp")
    create_ext_grid(net, bus=busn, vm_pu=1.0, name="Grid Connection", s_sc_max_mva=5000,
                    rx_max=0.1, r0x0_max=0.1, x0x_max=1.0)
    create_std_type(net, {"r0_ohm_per_km": .154, "x0_ohm_per_km": 0.5277876,
                          "c0_nf_per_km": 170.4, "g0_us_per_km": 0, "max_i_ka": 0.741,
                          "r_ohm_per_km": .0385, "x_ohm_per_km": 0.1319469,
                          "c_nf_per_km": 170}, "example_type3")
    create_line(net, from_bus=busn, to_bus=busm, length_km=1.0, std_type="example_type3")
    create_line(net, from_bus=busn, to_bus=busp, length_km=1.0, std_type="example_type3")
    create_line_from_parameters(net, from_bus=busn, to_bus=busk, length_km=1.0, r0_ohm_per_km=.1005,
                                x0_ohm_per_km=0.4900884, c0_nf_per_km=200.5, max_i_ka=0.89,
                                r_ohm_per_km=.0251, x_ohm_per_km=0.1225221, c_nf_per_km=210)
    create_line_from_parameters(net, from_bus=busk, to_bus=busm, length_km=1.0,
                                r0_ohm_per_km=0.0848, x0_ohm_per_km=0.4649556, c0_nf_per_km=230.6,
                                max_i_ka=0.963, r_ohm_per_km=0.0212, x_ohm_per_km=0.1162389, c_nf_per_km=230)
    create_line_from_parameters(net, from_bus=busk, to_bus=busp, length_km=1.0, r0_ohm_per_km=.3048,
                                x0_ohm_per_km=0.6031856, c0_nf_per_km=140.3, max_i_ka=0.531,
                                r_ohm_per_km=.0762, x_ohm_per_km=0.1507964, c_nf_per_km=140)
    add_zero_impedance_parameters(net)

    create_asymmetric_load(net, busk, p_a_mw=50, q_a_mvar=20, p_b_mw=80, q_b_mvar=60,
                           p_c_mw=20, q_c_mvar=5)
    create_asymmetric_load(net, busm, p_a_mw=50, q_a_mvar=50, p_b_mw=10, q_b_mvar=15,
                           p_c_mw=10, q_c_mvar=5)
    create_asymmetric_load(net, busp, p_a_mw=50, q_a_mvar=20, p_b_mw=60, q_b_mvar=20,
                           p_c_mw=10, q_c_mvar=5)
    runpp_3ph_with_consistency_checks(net, init=init, recycle=recycle)
    runpp_3ph_with_consistency_checks(net, init=init, recycle=recycle)
    assert net['converged']

    bus_pp = np.abs(net.res_bus_3ph[['vm_a_pu', 'vm_b_pu', 'vm_c_pu']]
                    [~np.isnan(net.res_bus_3ph.vm_a_pu)].values)
    bus_pf = np.abs(np.array([[0.98085729, 0.97711997, 1.04353786],
                              [0.97828577, 0.97534651, 1.04470864],
                              [0.97774307, 0.97648197, 1.04421233],
                              [0.9780892, 0.97586805, 1.04471106]]))
    assert np.max(np.abs(bus_pp - bus_pf)) < 1e-8

    line_pp = np.abs(net.res_line_3ph[
                         ['i_a_from_ka', 'i_b_from_ka', 'i_c_from_ka',
                          'i_a_to_ka', 'i_b_to_ka', 'i_c_to_ka',
                          'p_a_from_mw', 'p_b_from_mw', 'p_c_from_mw',
                          'q_a_from_mvar', 'q_b_from_mvar', 'q_c_from_mvar',
                          'p_a_to_mw', 'p_b_to_mw', 'p_c_to_mw',
                          'q_a_to_mvar', 'q_b_to_mvar', 'q_c_to_mvar',
                          'loading_a_percent', 'loading_b_percent', 'loading_c_percent',
                          'loading_percent']].values)
    line_pf = np.abs(np.array(
        [[0.98898804851, 0.68943734, 0.19848961,
          0.99093993, 0.69146384, 0.19966503,
          49.87434308, 33.86579548, 12.44659879,
          36.16562613, 26.14426519, 4.25746428,
          -49.75842138, -33.90236497, -12.45155362,
          -36.19862688, -26.25675246, -4.50384238,
          133.730100000000, 93.314960000000, 26.945350000000,
          133.730100000000],
         [0.87075816277, 1.03463205, 0.19072622,
          0.87210779, 1.03599167, 0.19188991,
          49.59359423, 58.53676842, 11.97553941,
          21.96967200, 26.37559958, 4.04458873,
          -49.47110289, -58.55284705, -11.98669516,
          -22.07474008, -26.34476811, -4.29078447,
          117.693400000000, 139.809900000000, 25.896070000000,
          139.809900000000],
         [0.95760407055, 1.14786582, 0.24829126,
          0.95975383, 1.15028040, 0.24975553,
          50.87938854, 57.53628873, 15.54470531,
          31.13888557, 41.99378843, 5.39758513,
          -50.76249094, -57.56374777, -15.56099267,
          -31.28560646, -41.99056453, -5.69609575,
          107.837500000000, 129.245000000000, 28.062420000000,
          129.245000000000],
         [0.21780921494, 0.42795803, 0.03706412,
          0.22229619, 0.42603286, 0.03771703,
          0.23292404, -23.88471674, -2.45255095,
          13.53037092, -11.49972060, 0.17971665,
          -0.24157862, 23.90236497, 2.45155361,
          -13.80137312, 11.25675247, -0.49615762,
          23.083720000000, 44.440090000000, 3.916618000000,
          44.440090000000],
         [0.03712221482, 0.10766244, 0.03093505,
          0.03446871, 0.10500386, 0.03179428,
          0.52956690, 1.44846452, -1.98645639,
          -2.24476446, -6.50971485, 0.51637910,
          -0.52889712, -1.44715295, 1.98669515,
          2.07474008, 6.34476812, -0.70921554,
          6.991001000000, 20.275410000000, 5.987624000000,
          20.275410000000]]))
    assert np.max(np.abs(line_pp - line_pf)) < 1e-4


def test_3ph_bus_mapping_order():
    net = create_empty_network()
    b2 = create_bus(net, vn_kv=0.4, index=4)
    create_bus(net, vn_kv=0.4, in_service=False, index=3)
    b1 = create_bus(net, vn_kv=0.4, index=7)

    create_ext_grid(net, b1, vm_pu=1.0, s_sc_max_mva=10, rx_max=0.1)
    net.ext_grid["x0x_max"] = 1.
    net.ext_grid["r0x0_max"] = 0.1
    create_std_type(net, {"r_ohm_per_km": 0.1013, "x_ohm_per_km": 0.06911504,
                          "c_nf_per_km": 690, "g_us_per_km": 0, "max_i_ka": 0.44,
                          "c0_nf_per_km": 312.4, "g0_us_per_km": 0, "r0_ohm_per_km": 0.4053,
                          "x0_ohm_per_km": 0.2764602}, "N2XRY 3x185sm 0.6/1kV")

    create_line(net, b1, b2, 1.0, std_type="N2XRY 3x185sm 0.6/1kV", index=4)
    create_line(net, b1, b2, 1.0, std_type="N2XRY 3x185sm 0.6/1kV", index=3, in_service=False)
    create_line(net, b1, b2, 1.0, std_type="N2XRY 3x185sm 0.6/1kV", index=7)
    add_zero_impedance_parameters(net)
    create_load(net, b2, p_mw=0.030, q_mvar=0.030)
    runpp(net)
    runpp_3ph_with_consistency_checks(net)
    assert net['converged']

    assert np.allclose(net.res_bus_3ph.vm_a_pu.values, net.res_bus.vm_pu.values, equal_nan=True)
    assert net.res_bus_3ph.index.tolist() == net.res_bus.index.tolist()

    assert net.res_line_3ph.index.tolist() == net.res_line.index.tolist()
    assert np.allclose(net.res_line.p_from_mw, net.res_line_3ph.p_a_from_mw +
                       net.res_line_3ph.p_b_from_mw +
                       net.res_line_3ph.p_c_from_mw)
    assert np.allclose(net.res_line.loading_percent, net.res_line_3ph.loading_a_percent)


def test_3ph_two_bus_line_powerfactory():
    net = create_empty_network()

    b1 = create_bus(net, vn_kv=0.4)
    b2 = create_bus(net, vn_kv=0.4)

    create_ext_grid(net, b1, vm_pu=1.0, s_sc_max_mva=10, rx_max=0.1)
    net.ext_grid["x0x_max"] = 1.
    net.ext_grid["r0x0_max"] = 0.1
    create_std_type(net, {"r_ohm_per_km": 0.1013, "x_ohm_per_km": 0.06911504,
                          "c_nf_per_km": 690, "g_us_per_km": 0, "max_i_ka": 0.44,
                          "c0_nf_per_km": 312.4, "g0_us_per_km": 0, "r0_ohm_per_km": 0.4053,
                          "x0_ohm_per_km": 0.2764602}, "N2XRY 3x185sm 0.6/1kV")

    create_line(net, b1, b2, 0.4, std_type="N2XRY 3x185sm 0.6/1kV")
    add_zero_impedance_parameters(net)
    create_load(net, b2, p_mw=0.010, q_mvar=0.010)
    create_asymmetric_load(net, b2, p_a_mw=0.020, q_a_mvar=0.010, p_b_mw=0.015, q_b_mvar=0.005, p_c_mw=0.025,
                           q_c_mvar=0.010)
    runpp_3ph_with_consistency_checks(net)
    assert net['converged']

    bus_pp = np.abs(net.res_bus_3ph[['vm_a_pu', 'vm_b_pu', 'vm_c_pu']].values)
    bus_pf = np.abs(np.array([[0.99939853552, 1.0013885141, 0.99921580141],
                              [0.97401782343, 0.98945593737, 0.96329605983]]))

    assert np.max(np.abs(bus_pp - bus_pf)) < 4e-6

    line_pp = np.abs(net.res_line_3ph[
                         ['i_a_from_ka', 'i_b_from_ka', 'i_c_from_ka',
                          'i_a_to_ka', 'i_b_to_ka', 'i_c_to_ka',
                          'p_a_from_mw', 'p_b_from_mw', 'p_c_from_mw',
                          'q_a_from_mvar', 'q_b_from_mvar', 'q_c_from_mvar',
                          'p_a_to_mw', 'p_b_to_mw', 'p_c_to_mw',
                          'q_a_to_mvar', 'q_b_to_mvar', 'q_c_to_mvar']].values)
    line_pf = np.abs(np.array(
        [[0.11946088987, 0.08812337783, 0.14074226065,
          0.1194708224, 0.088131567331, 0.14075063601,
          0.023810539354, 0.01855791658, 0.029375192747,
          0.013901720672, 0.008421814704, 0.013852398586,
          -0.023333142958, -0.018333405987, -0.028331643666,
          -0.013332756527, -0.008333413919, -0.013332422725]]))
    assert np.max(np.abs(line_pp - line_pf)) < 1e-5

    line_load_pp = np.abs(net.res_line_3ph[
                              ['loading_a_percent', 'loading_b_percent', 'loading_c_percent',
                               'loading_percent']].values)
    line_load_pf = np.abs(np.array(
        [[27.1525, 20.0299, 31.98878,
          31.98878]]))
    assert np.max(np.abs(line_load_pp - line_load_pf)) < 1e-2


def check_bus_voltages(net, result, trafo_vector_group):
    res_vm_pu = []
    ordered_bus_table = net.bus.sort_values(by='name').reset_index(drop=True)
    for i in range(len(ordered_bus_table)):
        index = net.bus[net.bus.name == ordered_bus_table['name'][i]].index[0]
        res_vm_pu.append(net.res_bus_3ph.vm_a_pu[index])
        res_vm_pu.append(net.res_bus_3ph.vm_b_pu[index])
        res_vm_pu.append(net.res_bus_3ph.vm_c_pu[index])

    # max_tol = 0
    # for tol in [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7]:
    #     if np.allclose(result, res_vm_pu, atol=tol):
    #         max_tol = tol
    # print('Voltage Magnitude for %s is within tolerance of %s' % (trafo_vector_group, max_tol))

    tolerances = {'YNyn': 1e-05,
                  'Dyn': 1e-05,
                  'Yzn': 1e-03}

    assert np.allclose(result, res_vm_pu, atol=tolerances[trafo_vector_group],
                       rtol=0), f"Incorrect results for {trafo_vector_group}"


def check_line_currents(net, result, trafo_vector_group):
    res_line_i_ka = []
    ordered_line_table = net.line.sort_values(by='name').reset_index(drop=True)
    for i in range(len(ordered_line_table)):
        index = net.line[net.line.name == ordered_line_table['name'][i]].index[0]
        res_line_i_ka.append(net.res_line_3ph.i_a_from_ka[index])
        res_line_i_ka.append(net.res_line_3ph.i_b_from_ka[index])
        res_line_i_ka.append(net.res_line_3ph.i_c_from_ka[index])
        res_line_i_ka.append(net.res_line_3ph.i_a_to_ka[index])
        res_line_i_ka.append(net.res_line_3ph.i_b_to_ka[index])
        res_line_i_ka.append(net.res_line_3ph.i_c_to_ka[index])

    # max_tol = 0
    # for tol in [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7]:
    #     if np.allclose(result, res_line_i_ka, atol=tol):
    #         max_tol = tol
    # print('Line current for %s is within tolerance of %s' % (trafo_vector_group, max_tol))

    tolerances = {'YNyn': 1e-07,
                  'Dyn': 1e-07,
                  'Yzn': 1e-04}

    assert np.allclose(result, res_line_i_ka, atol=tolerances[trafo_vector_group])
    if not np.allclose(result, res_line_i_ka, atol=tolerances[trafo_vector_group]):
        raise ValueError("Incorrect results for vector group %s" % trafo_vector_group, res_line_i_ka, result)


def check_trafo_currents(net, result, trafo_vector_group):
    res_trafo_i_ka = []
    ordered_trafo_table = net.trafo.sort_values(by='name').reset_index(drop=True)
    for i in range(len(ordered_trafo_table)):
        index = net.trafo[net.trafo.name == ordered_trafo_table['name'][i]].index[0]
        res_trafo_i_ka.append(net.res_trafo_3ph.i_a_hv_ka[index])
        res_trafo_i_ka.append(net.res_trafo_3ph.i_b_hv_ka[index])
        res_trafo_i_ka.append(net.res_trafo_3ph.i_c_hv_ka[index])
        res_trafo_i_ka.append(net.res_trafo_3ph.i_a_lv_ka[index])
        res_trafo_i_ka.append(net.res_trafo_3ph.i_b_lv_ka[index])
        res_trafo_i_ka.append(net.res_trafo_3ph.i_c_lv_ka[index])

    # max_tol = 0
    # for tol in [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7]:
    #     if np.allclose(result, res_trafo_i_ka, atol=tol):
    #         max_tol = tol
    # print('Trafo current for %s is within tolerance of %s' % (trafo_vector_group, max_tol))

    tolerances = {'YNyn': 1e-03,
                  'Dyn': 1e-03,
                  'Yzn': 1e-03}

    assert np.allclose(result, res_trafo_i_ka, atol=tolerances[trafo_vector_group])
    if not np.allclose(result, res_trafo_i_ka, atol=tolerances[trafo_vector_group]):
        raise ValueError("Incorrect results for vector group %s" % trafo_vector_group, res_trafo_i_ka, result)


def check_results(net, trafo_vector_group, results):
    check_bus_voltages(net, results[0], trafo_vector_group)
    check_line_currents(net, results[1], trafo_vector_group)
    check_trafo_currents(net, results[2], trafo_vector_group)


def make_nw(net, bushv, tap_ps, case, vector_group):
    b1 = create_bus(net, bushv, zone=vector_group, index=get_free_id(net.bus))
    b2 = create_bus(net, 0.4, zone=vector_group)
    b3 = create_bus(net, 0.4, zone=vector_group)
    create_ext_grid(net, b1, s_sc_max_mva=10000,
                    rx_max=0.1, r0x0_max=0.1, x0x_max=1.0)
    create_transformer_from_parameters(net, hv_bus=b1, lv_bus=b2,
                                       sn_mva=1.6, vn_hv_kv=10,
                                       vn_lv_kv=0.4, vk_percent=6,
                                       vkr_percent=0.78125, pfe_kw=2.7,
                                       i0_percent=0.16875, shift_degree=0,
                                       tap_side='lv', tap_neutral=0,
                                       tap_min=-2, tap_max=2,
                                       tap_step_degree=0,
                                       tap_step_percent=2.5,
                                       tap_changer_type="Ratio",
                                       vk0_percent=6, vkr0_percent=0.78125,
                                       mag0_percent=100, mag0_rx=0.,
                                       si0_hv_partial=0.9, vector_group=vector_group,
                                       parallel=1, tap_pos=tap_ps,
                                       index=get_free_id(net.trafo) + 1)
    create_line_from_parameters(net, b2, b3, length_km=0.5, r_ohm_per_km=0.1941, x_ohm_per_km=0.07476991,
                                c_nf_per_km=1160., max_i_ka=0.421,
                                endtemp_degree=70.0, r0_ohm_per_km=0.7766,
                                x0_ohm_per_km=0.2990796,
                                c0_nf_per_km=496.2,
                                index=get_free_id(net.line) + 1)
    if case == "bal_wye":
        # Symmetric Load
        create_load(net, b3, 0.08, 0.012, type='wye')
    elif case == "delta_wye":
        # Unsymmetric Light Load
        create_asymmetric_load(net, b3, p_a_mw=0.0044, q_a_mvar=0.0013, p_b_mw=0.0044, q_b_mvar=0.0013,
                               p_c_mw=0.0032, q_c_mvar=0.0013, type='wye')
        create_asymmetric_load(net, b3, p_a_mw=0.0300, q_a_mvar=0.0048, p_b_mw=0.0280, q_b_mvar=0.0036,
                               p_c_mw=0.027, q_c_mvar=0.0043, type='delta')
    elif case == "wye" or case == "delta":
        # Unsymmetric Heavy Load
        create_asymmetric_load(net, b3, p_a_mw=0.0300, q_a_mvar=0.0048, p_b_mw=0.0280, q_b_mvar=0.0036,
                               p_c_mw=0.027, q_c_mvar=0.0043, type=case)


def test_trafo_asym():
    nw_dir = os.path.abspath(os.path.join(pp_dir, "test/loadflow"))
    # only 3 vector groups are supported in the 3ph power flow
    for trafo_vector_group in ["YNyn", "Dyn", "Yzn"]:
        net = from_json(nw_dir + '/runpp_3ph Validation.json')
        net['trafo'].vector_group = trafo_vector_group
        runpp_3ph_with_consistency_checks(net)
        assert net['converged']
        check_results(net, trafo_vector_group, get_PF_Results(trafo_vector_group))

def _test_trafo_shifts(net, rtol):
    # Dyn
    for clock in [-30, 30, 150, 210, -150]:
        net.trafo.vector_group = "Dyn"
        net.trafo.shift_degree = clock
        runpp_3ph_with_consistency_checks(net)
        trafo_currents_consistent_3ph(net, rtol)

    # YNyn
    for clock in [0, 180, -180]:
        net['trafo'].vector_group = "YNyn"
        net['trafo'].shift_degree = clock
        runpp_3ph_with_consistency_checks(net)
        trafo_currents_consistent_3ph(net, rtol)

def test_trafo_asym_currents__high_neg_seq():
    """
    Tests trafo currents consistency in asymmetric load flow.
    In order to verify current consistency between HT and LT side, we need to make
    shunt admittance zero, so that no load currents do not disturb HT side currents.
    For this reason, pfe_kw = 0, i0_percent = 0, mag0_percent = BIG_NUMBER.
    """
    net = create_empty_network()
    add_zero_impedance_parameters(net)
    create_bus(net, 11, "source")
    create_bus(net, 11, "HT")
    create_bus(net, 0.4, "LT")
    create_ext_grid(net, bus=0, vm_pu=1.0, s_sc_max_mva=9.99e20, rx_max=0.1, x0x_max=1.0, r0x0_max=0.1)
    create_line_from_parameters(net, 0, 1, 1, 0.2, 0.33, 0, max_i_ka=10, r0_ohm_per_km=0.2, x0_ohm_per_km=0.33,
                                   c0_nf_per_km=0, g0_gs_per_km=0)
    create_transformer_from_parameters(net=net, hv_bus=1, lv_bus=2, sn_mva=10, vn_hv_kv=11, vn_lv_kv=0.415,
                                          vkr_percent=0, vk_percent=4, pfe_kw=0, i0_percent=0,
                                          shift_degree=-30, vector_group="Dyn",
                                          vk0_percent=4, vkr0_percent=0, mag0_percent=1e20, mag0_rx=0, si0_hv_partial=50)
    create_asymmetric_load(net, 2, 1.6, 0.44019238, 0.95980762, 0.8, 0.86961524, -0.16961524)
    _test_trafo_shifts(net, rtol=1e-9)

def test_2trafos():
    net = create_empty_network()
    make_nw(net, 10., 0., "wye", "YNyn")
    make_nw(net, 10., 0., "wye", "YNyn")
    runpp_3ph_with_consistency_checks(net)
    assert net['converged']
    assert np.allclose(net.res_ext_grid_3ph.iloc[0].values, net.res_ext_grid_3ph.iloc[1].values)


def test_3ph_isolated_nodes():
    v_base = 110  # 110kV Base Voltage
    mva_base = 100  # 100 MVA
    net = create_empty_network(sn_mva=mva_base)

    busn = create_bus(net, vn_kv=v_base, name="busn", index=1)
    create_bus(net, vn_kv=20., in_service=True, index=2, name="busx")
    busk = create_bus(net, vn_kv=v_base, name="busk", index=5)
    busl = create_bus(net, vn_kv=v_base, name="busl", index=6)
    create_bus(net, vn_kv=20., in_service=False, index=3)
    busy = create_bus(net, vn_kv=20., in_service=True, index=0, name="busy")

    create_ext_grid(net, bus=busn, vm_pu=1.0, name="Grid Connection",
                    s_sc_max_mva=5000, rx_max=0.1)
    net.ext_grid["r0x0_max"] = 0.1
    net.ext_grid["x0x_max"] = 1.0
    create_std_type(net, {"r0_ohm_per_km": 0.0848, "x0_ohm_per_km": 0.4649556,
                          "c0_nf_per_km": 230.6, "g0_us_per_km": 0, "max_i_ka": 0.963,
                          "r_ohm_per_km": 0.0212, "x_ohm_per_km": 0.1162389,
                          "c_nf_per_km": 230}, "example_type")
    # Loads on supplied buses
    create_asymmetric_load(net, busk, p_a_mw=50, q_a_mvar=50, p_b_mw=10, q_b_mvar=15,
                           p_c_mw=10, q_c_mvar=5)
    create_load(net, bus=busl, p_mw=7, q_mvar=0.070, name="Load 1")
    # Loads on unsupplied buses
    create_load(net, bus=busy, p_mw=70, q_mvar=70, name="Load Y")
    create_line(net, from_bus=busn, to_bus=busk, length_km=50.0, std_type="example_type")
    create_line(net, from_bus=busl, to_bus=busk, length_km=50.0, std_type="example_type")
    add_zero_impedance_parameters(net)
    runpp_3ph_with_consistency_checks(net)
    assert net['converged']
    assert np.allclose(net.res_bus_3ph.T[[0, 2, 3]].T[["vm_a_pu", "va_a_degree", "vm_b_pu",
                                                       "va_b_degree", "vm_c_pu", "va_c_degree"]], np.nan,
                       equal_nan=True)
    assert np.allclose(net.res_bus_3ph.T[[0, 2, 3]].T[["p_a_mw", "q_a_mvar", "p_b_mw", "q_b_mvar",
                                                       "p_c_mw", "q_c_mvar"]], 0.0)


def test_balanced_power_flow_with_unbalanced_loads_and_sgens():
    net = create_empty_network(sn_mva=100)
    make_nw(net, 10, 0, "wye", "Dyn")
    create_asymmetric_sgen(net, 1, p_a_mw=0.01, p_b_mw=0.02, scaling=0.8)
    runpp_with_consistency_checks(net)

    vm_pu = net.res_bus.vm_pu

    net.asymmetric_load.in_service = False
    create_load(net, bus=net.asymmetric_load.bus.iloc[0],
                scaling=net.asymmetric_load.scaling.iloc[0],
                p_mw=net.asymmetric_load.loc[0, ["p_a_mw", "p_b_mw", "p_c_mw"]].sum(),
                q_mvar=net.asymmetric_load.loc[0, ["q_a_mvar", "q_b_mvar", "q_c_mvar"]].sum()
                )
    runpp_with_consistency_checks(net)
    assert net.res_bus.vm_pu.equals(vm_pu)

    net.asymmetric_sgen.in_service = False
    create_sgen(net, bus=net.asymmetric_sgen.bus.iloc[0],
                scaling=net.asymmetric_sgen.scaling.iloc[0],
                p_mw=net.asymmetric_sgen.loc[0, ["p_a_mw", "p_b_mw", "p_c_mw"]].sum(),
                q_mvar=net.asymmetric_sgen.loc[0, ["q_a_mvar", "q_b_mvar", "q_c_mvar"]].sum()
                )
    runpp_with_consistency_checks(net)
    assert net.res_bus.vm_pu.equals(vm_pu)


def test_3ph_with_impedance():
    nw_dir = os.path.abspath(os.path.join(pp_dir, "test/loadflow"))
    net = from_json(nw_dir + '/runpp_3ph Validation.json')
    net.line.c_nf_per_km = 0.
    net.line.g_nf_per_km = 0.
    net.line["c0_nf_per_km"] = 0.
    net.line["g0_us_per_km"] = 0.
    net_imp = copy.deepcopy(net)
    replace_line_by_impedance(net_imp, net.line.index, 100)
    runpp_3ph(net)
    runpp_3ph(net_imp)
    assert dataframes_equal(net.res_bus_3ph, net_imp.res_bus_3ph)


def test_shunt_3ph():
    net = create_empty_network()
    b1 = create_bus(net, 20.0)
    b2 = create_bus(net, 20.0)
    create_ext_grid(net, b1, s_sc_max_mva=1000, rx_max=0.1, x0x_max=1.0, r0x0_max=0.1)
    create_line_from_parameters(net, b1, b2, length_km=1, r_ohm_per_km=1, x_ohm_per_km=1, c_nf_per_km=1,
                                r0_ohm_per_km=1, x0_ohm_per_km=1, c0_nf_per_km=1, max_i_ka=1)
    create_shunt(net, in_service=True, bus=b2, p_mw=1, q_mvar=1)
    runpp_3ph(net)

    bus_cols = ["vm_a_pu", "va_a_degree", "vm_b_pu", "va_b_degree", "vm_c_pu", "va_c_degree"]
    bus_pp = np.abs(net.res_bus_3ph[bus_cols].values)
    bus_expected = np.abs(np.array([[1.0, 0.0, 1.0, -120.0, 1.0, 120.0],
                                    [0.9950250311521571, -8.955222949441059e-06, 0.995025031152157, -120.00000895522294,
                                     0.995025031152157, 119.99999104477705]]))

    assert np.max(np.abs(bus_pp - bus_expected)) < 1e-5

    line_cols = ["p_a_from_mw", "q_a_from_mvar", "p_b_from_mw", "q_b_from_mvar", "q_c_from_mvar", "p_a_to_mw",
                 "q_a_to_mvar", "p_b_to_mw", "q_b_to_mvar", "p_c_to_mw", "q_c_to_mvar"]
    line_power_pp = np.abs(net.res_line_3ph[line_cols].values)
    line_power_expected = np.abs(np.array(
        [[0.33167495789351165, 0.33163327786948577, 0.33167495789351165, 0.33163327786948565, 0.33163327786948565,
          -0.3300249368894127, -0.3300249368948133, -0.3300249368894127, -0.3300249368948133, -0.3300249368894127,
          -0.3300249368948134]]))
    assert np.max(np.abs(line_power_pp - line_power_expected)) < 1e-5


@pytest.mark.xfail
def test_3ph_enforce_q_lims(result_test_network):
    v_tol = 1e-6
    s_tol = 5e-3

    net = result_test_network
    buses = net.bus[net.bus.zone == "test_enforce_qlims"]
    gens = [x for x in net.gen.index if net.gen.bus[x] in buses.index]
    b2 = buses.index[1]
    b3 = buses.index[2]
    g1 = gens[0]

    # enforce reactive power limits
    runpp_3ph(net, enforce_q_lims=True)

    # powerfactory results
    u2 = 1.00607194
    u3 = 1.00045091

    assert abs(net.res_bus.vm_pu.at[b2] - u2) < v_tol
    assert abs(net.res_bus.vm_pu.at[b3] - u3) < v_tol
    assert abs(net.res_gen.q_mvar.at[g1] - net.gen.min_q_mvar.at[g1]) < s_tol


@pytest.mark.xfail
def test_recycle_pq():
    net = create_empty_network()
    _, b2, _ = add_grid_connection(net)
    pl = 1.2
    ql = 1.1
    ps = 0.5
    u_set = 1.0

    b3 = create_bus(net, vn_kv=0.4)
    create_bus(net, vn_kv=0.4, in_service=False)
    create_line_from_parameters(
        net, b2, b3, 12.2, r_ohm_per_km=0.08, x_ohm_per_km=0.12, c_nf_per_km=300, max_i_ka=0.2, df=0.8
    )
    create_load(net, b3, p_mw=pl, q_mvar=ql)
    create_gen(net, b2, p_mw=ps, vm_pu=u_set)
    pl = 1.2
    ql = 0.0
    net["load"].at[0, "q_mvar"] = ql
    runpp_3ph(net, recycle={"trafo": False, "gen": False, "bus_pq": True})
    assert np.allclose(net.res_load.at[0, "p_mw"], pl)
    assert np.allclose(net.res_load.at[0, "q_mvar"], ql)
    pl = 0.8
    ql = 0.55
    net.load.at[0, "p_mw"] = pl
    net["load"].at[0, "q_mvar"] = ql
    runpp_3ph(net, recycle={"bus_pq": True, "trafo": False, "gen": False})
    assert np.allclose(net.res_load.p_mw.iloc[0], pl)
    assert np.allclose(net.res_load.q_mvar.iloc[0], ql)


@pytest.mark.xfail
def test_connectivity_check_island_without_pv_bus():
    # Network with islands without pv bus -> all buses in island should be set out of service
    net = create_cigre_network_mv(with_der=False)
    iso_buses, iso_p, iso_q, *_ = get_isolated(net)
    assert len(iso_buses) == 0
    assert np.isclose(iso_p, 0)
    assert np.isclose(iso_q, 0)

    isolated_bus1 = create_bus(net, vn_kv=20.0, name="isolated Bus1")
    isolated_bus2 = create_bus(net, vn_kv=20.0, name="isolated Bus2")
    create_line(
        net, isolated_bus2, isolated_bus1, length_km=1, std_type="N2XS(FL)2Y 1x300 RM/35 64/110 kV", name="IsolatedLine"
    )
    iso_buses, iso_p, iso_q, *_ = get_isolated(net)
    assert len(iso_buses) == 2
    assert np.isclose(iso_p, 0)
    assert np.isclose(iso_q, 0)

    create_load(net, isolated_bus1, p_mw=0.2, q_mvar=0.02)
    create_sgen(net, isolated_bus2, p_mw=0.15, q_mvar=0.01)

    # with pytest.warns(UserWarning):
    iso_buses, iso_p, iso_q, *_ = get_isolated(net)
    assert len(iso_buses) == 2
    assert np.isclose(iso_p, 350)
    assert np.isclose(iso_q, 30)
    # with pytest.warns(UserWarning):
    runpp_3ph(net, check_connectivity=True)


@pytest.mark.xfail
@pytest.mark.parametrize("numba", [True, False])
def test_z_switch(numba):
    net = create_empty_network()
    for i in range(3):
        create_bus(net, vn_kv=0.4)
        create_load(net, i, p_mw=0.1)
    create_ext_grid(net, 0, vm_pu=1.0)
    create_line_from_parameters(
        net, 0, 1, 1, r_ohm_per_km=0.1 / np.sqrt(2), x_ohm_per_km=0.1 / np.sqrt(2), c_nf_per_km=0, max_i_ka=0.2
    )
    create_switch(net, 0, 2, et="b", z_ohm=0.1)

    runpp_3ph(net, numba=numba, switch_rx_ratio=1)
    assert pytest.approx(net.res_bus.vm_pu.at[1], abs=1e-9) == net.res_bus.vm_pu.at[2]

    net_zero_z_switch = copy.deepcopy(net)
    net_zero_z_switch.switch.z_ohm = 0
    runpp(net_zero_z_switch, numba=numba, switch_rx_ratio=1)
    assert pytest.approx(net_zero_z_switch.res_bus.vm_pu.at[0], abs=1e-9) == net_zero_z_switch.res_bus.vm_pu.at[2]


if __name__ == "__main__":
    pytest.main([__file__, "-xs"])
